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$-H ' Abstract. The phase diagram of a biased graphene bilayer is computed and the 

(/3 , existence of a ferromagnetic phase is discussed both in the critical on-site interaction 

Uc versus doping density and versus temperature. We show that in the ferromagnetic 
phase the two planes have unequal magnetization and that the electronic density is 
hole like in one plane and electron like in the other. We give evidence for a first-order 
' O ' phase transition between paramagnetic and ferromagnetic phases induced by doping 

^ , at zero temperature. 
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First order ferromagnetic phase transition in the low electronic density regime of a biased graphene bilayer2 
1. Introduction 

Graphene, a two-dimensional hexagonal lattice of carbon atoms, has attracted 
considerable attention due to its unusual electronic properties, characterized by massless 
Dirac Fermions.fl], O [3j It was first produced via micromechanical cleavage on top of a 
Si02 substrate[H [5| and its hallmark is the half integer quantum Hall effect. [HIT] 

In addition to graphene, few-layer graphene can also be produced. Of particular 
interest to us is the double layer graphene system, where one encounters two carbon 
layers placed on top of each other according to usual Bernal stacking of graphite 
(see Fig. [1]). The low-energy properties of this so-called bilayer graphene are then 
described by massive Dirac Fermions.[8j These new quasi-particles have a quadratic 
dispersion close to the neutrality point and have recently been identified in Quantum 
Hall measurements [9j and in Raman spectroscopy. [TOl fTT] 

In a graphene bilayer it is possible to have the two planes at different electrostatic 
potentials. [Ill [I3] As a consequence, a gap opens up at the Dirac point and the low 
energy band acquires a Mexican hat relation dispersion. [T4] This system is called a 
biased graphene bilayer. The potential difference created between the two layers can 
be obtained by applying a back gate voltage to the bilayer system and covering the 
exposed surface with some chemical dopant, as for example Potassium [T2] or NHq.[T3] 
In addition, it is also possible to control the potential difference between the layers 
by using back and top gate setups. [T5] The opening of the gap at the Dirac point 
in the biased bilayer system was demonstrated both by angle resolved photoemission 
experiments (ARPES)|l2| and Hall effect measurements. [T3] The electronic gap in the 
biased system has been also observed in epitaxially grown graphene films on SiC crystal 
surfaces. 1 16) 

Due to the Mexican hat dispersion relation the density of states close to the gap 
diverges as the square root of the energy. The possibility of having an arbitrary large 
density of states at the Fermi energy poses the question whether this system can be 
unstable toward a ferromagnetic ground state. The question of magnetism in carbon 
based systems already has a long history. Even before the discovery of graphene, 
highly oriented pyrolytic graphite (HOPG) has attracted a broad interest due to the 
observation of anomalous properties, such as magnetism and insulating behavior in the 
direction perpendicular to the planes. [ITl [181 [El [201 [21] The research oi s — p based 
magnetism[22|, [23], [24] was especially motivated by the technological use of nanosized 
particles of graphite, which show interesting features depending on their shape, edges, 
and applied field of pressure. [25] 

Microscopic theoretical models of bulk carbon magnetism include nitrogen-carbon 
compositions where ferromagnetic ordering of spins could exist in tt delocalized systems 
due to a lone electron pair on a trivalent element [26] or intermediate graphite- 
diamond structures where the alternating sp"^ and sp^ carbon atoms play the role of 
different valence elements. [27] More general models focus on the interplay between 
disorder and interaction. ^8], [29] Further, midgap states due to zig-zag edges play a 
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predominant role in the formation of magnetic moments [301 [31] which support flat- 
band ferromagnetism.[32 l [33l [34] A generic model based on midgap states was recently 
proposed in Refs. [35|, [36]. Magnetism is also found in fullerene based metal- free 
systems. [37] For a recent overview on metal- free Carbon-based magnetism see Ref. [38) . 

To understand carbon-based magnetism in graphite, one may start with the 
simplest case of one-layer, i.e., graphene. Because the density of states of intrinsic 
graphene vanishes at the Dirac point, the simple Stoner-like argument predicts an 
arbitrary large value of the Coulomb on-site energy needed to produce a ferromagnetic 
ground state. [39l [40] In fact, because of the vanishing density of states, the Coulomb 
interaction is not screened and the Hubbard model is not a good starting point to 
study ferromagnetism in clean graphene. One, therefore, has to consider the exchange 
instability of the Dirac gas due to the bare, long-range Coulomb interaction in two 
dimensions. This study shows that for a clean, doped or undoped graphene layer, a spin- 
polarized ground-state due to the gain of exchange energy is only favorable for unphysical 
values of the dimensionless coupling constant of graphene. [41] The paramagnetic ground- 
state of clean graphene is thus stable against the exchange interaction. If the system is 
disordered, e.g., due to vacancies or edge states, a flnite density of states builds in at 
the Dirac point. As a consequence, a flnite Hubbard interaction for driving the system 
to a ferromagnetic ground state is obtained. [42] In this case, the exchange interaction 
favors a ferromagnetic ground state for reasonable values of the dimensionless coupling 
parameter. [H] The presence of itinerant magnetism due to quasi-localized states induced 
by single-atom defects in graphene, such as vacancies, [43] has also been obtained recently 
using flrst-principles. [44] 

The situation is quite different in a bilayer system. There, a flnite density of 
states exists at the neutrality point producing some amount of screening in the system. 
Moreover, in the case of a biased bilayer and for densities close to the energy gap, the 
density of states is very large producing very effective screening. As a consequence, for 
this system the Hubbard model is a good starting point to study the tendency toward 
ferromagnetism. From the point of view of the exchange instability of the bilayer system, 
it is found that the system is always unstable toward a ferromagnetic ground state for 
low enough particle densities. [45l [46l [47] 

In the present paper, we want to explore the fact that the Hubbard model is a 
good starting point to describe the Coulomb interactions in the regime where the Fermi 
energy is close to the band edge of the biased bilayer system. In particular we want to 
study the phase diagram of the system as function of the doping. We further want to 
determine the mean held critical temperature. 

The paper is organized as followed. In section [21 we introduce the model and deflne 
the mean-fleld decoupling which allows for different electronic density and magnetization 
in the two layers. In section [3l we set up the mean-fleld equations and present the 
numerical results in section [H We close with conclusions and future research directions. 
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2. Model Hamiltonian and mean field approximation 

The Hamiltonian of a biased bilayer Hubbard model is the sum of two pieces H = 
Htb + Hu^ where Htb is the tight-binding part and Hu is the Coulomb on-site 
interaction part of the Hamiltonian. The tight-binding Hamiltonian is itself a sum 
of four terms describing the tight-binding Hamiltonian of each plane, the hopping term 
between the planes, and the electrostatic bias applied to the two planes. We therefore 
have, 

2 

Htb = 2_^ Htb,i. + H± + Hy , (1) 



with 



Htb,. = - tYM\AR)b..{R) + al{R)b,^{R - a^) 

R,<j 

+ al{R)bUR-ci2) + H.c.], (2) 

H± = -t± 5^[aL(i?)62.(i?) + bl{R)a,,{R)] , (3) 



R,a 



and 



V v^ 

Hv = J 2^[naiaiR) + ribiaiR) - na2aiR) - nb2aiR)] ■ (4) 

R,a 

As regards the bias term in Eq. ([4]), we assume here that V can be externally controlled 
and is independent of the charge density in the system. This situation can be realized 
with a back and top gate setup. p!5] The on-site Coulomb part is given by, 

Hu = U ^ Kit (i?)n„i|(i?) +72611 (i?)nfeii(i2) 

R 

+ na2\ {R)na2i (R) + ^621 (^)^fe2i (R)] , (5) 

where rixLo-iR) = ^laiR)^!-a{R), with x = a,b, t = 1,2 and a =t, i- 

The problem defined by the Hamiltonian Htb + Hu can not be solved exactly 
and therefore we have to rest upon some approximation. Here we adopt a mean field 
approach, neglecting quantum fluctuations. Since we are interested in studying the 
existence of a ferromagnetic phase we have to introduce a broken symmetry ground 
state. There is however an important point to remark: since the two planes of the 
bilayer are at different electrostatic potentials one should expect that the electronic 
density and the magnetization will not be evenly distributed among the two layers. 
Therefore our broken symmetry ground state must take this aspect into account. As a 
consequence we propose the following broken symmetry ground state: 

/ ^N > n + An m + Am . . 

{n,i^{R)) = + cr , (6) 



and 



/„N> n — An m — Am , , 

(ri.2.(i?)) = 5 + o 3 , (7) 
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Figure 1. (Color online) The unit cell a of graphene bilayer in the Bernal stacking. 
The dashed hexagons are on top of the solid ones. The unit cell vectors have coordinates 
ai — a(3, ■\/3)/2 and a2 = a(3, — -\/3)/2. 



where n is the density per unit cell and m = n^ — n| is the spin polarization per unit 
cell. The quantities An and Am represent the difference in the electronic density and in 
the spin polarization between the two layers, respectively. [48] We note that m and Am 
are independent parameters, being in principle possible to have a ground state where 
m = but Am ^ 0. 

When transformed to momentum space the mean field Hamiltonian obtained from 
the above reads 

Hmf = Yl *fc,<7^fc,-^fc- 

k,cr 

KU 



[{n + An) — (m + Am) 



32 

-[{n — An)^ — (m — Am)^] 



32 
with ^I,, = [a\kaAka^(^\kaAka\ ^nd Hk,a givcu by 

/ Sa -t(f)k -t± \ 

-t(pl s,, 

p^ -t0fc 

-ti. -t(t)l p„ J 

with s, = f + (^ - a^^^) U, p, = -^ + (n^ - a^2^) U, and 0^ 
-^ _l_ gjfcai _|_ gjfca2 fjij^g energy eigenvalues are given by. 



(8) 



H^ 



k,cr 



\ 



(9) 



Ei'Uk^m^Am,) 



n 



a^]u 



(10) 



+ y^tl + V^ + At^\M' + J2x/ti + mtl + V^)\M^ 



where l,j = ± and V„ is given by 

V^ = V + UAn - aUAm , 



(11) 
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where we have introduced the definitions An = 4An and Am = 4Am. It is clear that 
as long as An and Am are finite the system has an eff'ective Va that differs from the 
bare value V. The momentum values are given by, 

(12) 






with mi,m2 = 0, 1, . . . A^ — 1, the number of unit cells given by iV^ = iV^, and 6i and 
62 given by, 



61 = — 1,V3 , b2 = —il,-V3). 
Sa Sa 

The Brillouin zone of the system is represented in Fig. [2l 



(13) 




Figure 2. (Color online) Brillouin zone of the bilayer problem. The Dirac point K 
has coordinates 27r(l, \/3/3)/(3a) and the M point has coordinates 27r(l,0)/(3a). 



3. Free energy and mean-field equations 



The free energy per unit cell, /, of Hamiltonian (E]) is given by. 



/= - 



Nr 



k,cr l,j=± 



,-(£^-^(fc)-/.)/(fcsT) 



(14) 



U 



~ —\n^ -m^ + {Anf - (Am)^] + /xn , 
16 

where /i is the chemical potential. 

Let us introduce the density of states per spin per unit cell p{E) defined as 



p{E) = j^Y.^{E-t\<P^\). 



(15) 
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The momentum integral in Eq. (ITSl) is over the Brillouin zone defined in Fig. [21 using 
the momentum definition (fT2l) . The integral can be performed leading to, 



p{E) 



2-Tr2 



t% 



\^(^^\ , t<E <3t. 



(16) 



y/Wt V 4^/* 7 



where F(x) is given by 

F{x) = {l + xr-^-^^, (17) 

and K(m) is defined as, 

K(m)= / ds[(l-a;')(l-mx2)]-i/2. (18) 

Using Eq. f lTSl ). the free energy in Eq. ( ITS! ) can be written as a one-dimensional 
integral, 

/ = 

a l,j=±-^ 

- — [n^ -771^ + (An)2 - (Am)^] +/in. (19) 

16 

The mean field equations are now obtained from the minimization of the free energy 

( fT9l) . The doping, Sn, relative to the situation where the system is at half filling is 

defined as, 

^^ = E E / dEp{E)f[E';^{E) - /.] - 4 , (20) 

where f{x) = (1 + e^^^'''^'^^)~^. The spin polarization per unit cell obeys the mean field 
equation, 

m = E E ^ / dEp{E)f[El;\E) - /i] . (21) 

a l,j=± -^ 

The difference in the electronic density between the two layers is obtained from, 

^^ = ^E E [dEpiE)f[El;^{E)-p^]v^^^{E), (22) 

where f^-'(-E') is given by 

„U^E) = 1^ f 1 + ^'^' ] . (23) 

and 



2tl + V^ + 4^2 + j2Jt\ + AE\tl + V^) . (24) 
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The difference in the magnetization between the two layers is obtained from 

1 v^ v^ 

2 



^^ = ^1^1^ ^JdEp{E)f[El;\E)-M'{E), (25) 



a l,j = ±l 

Let us now assume that the system supports a ferromagnetic ground state whose 
magnetization vanishes at some critical value Uc at zero temperature. Additionally we 
assume that Am = when m = 0, which will be shown to be the case in this system. 
The value of Uc is determined from expanding (|2T]1 to first order in m, leading to, 



^ Y, ldEp{E)S[E^^^{E,0,0)~f,] 



i,j=±i 



"-itj^,m-E-,e(E; 



l,j,k=±l 

= -jPbif^, Uc) = Ucpbif^, Uc) , (26) 

where Pb{p^, Uc) is the density of states per unit cell per spin for a biased bilayer at the 
energy jl = ^ — nUc/S and Pfe(/i, Uc) is the density of states per spin per lattice point. 
Although Eq. (I26l) looks like the usual Stoner criterion the fact that the bias V^ given in 
Eq. (ITTl) depend on U due to the difference in the electronic density An makes Eq. (l26l) a 
non-linear equation for Uc which must be solved self-consistently. For low doping 6n the 
product UcAh is a small number when compared to V and therefore it can be neglected 
in Eq. (ITT]) . In this case Eq. (I26l l reduces to the usual Stoner criterion: 

Uc ^ l/[p.(/i)] . (27) 

The quantities El in Eq. (J26l) are the roots of the delta function argument, 

E^;^iEl)-^^ = 0. (28) 

The quantity fij{El) is the derivative in order to the energy E of Eq. (1281 ) evaluated at 
the roots E^.. The roots El are given by 



Et = ^y V + v;2 + k2^4j,^^tl + V^) - t\V^ , (29) 

with /c = ±. Equation (l28ll cannot be solved for all bands: the existence of a solution 
is determined by /i. As a consequence we added the * symbol in the summation of 
Eq. (I26l ). which means that only bands for which Eq. (I28l) can be solved (two at the 
most) contribute to the summation. It also means that for the contributing bands only 
real roots in Eq. (I29ll are taken into account to the summation. The number of real 
roots in Eq. (l29l) depends on the particular band an /i through Eq. (l28l) . As the function 
fi,j{E) is given by 



it is clear that both roots are imaginary for jl in the range 

< u < , 31 
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which means that the system has an energy gap of value 

A, = /^^" . (32) 

We finally note that since we have assumed Am = 0, V^ does not eff'ectively depend on 
a. 

4. Results and Discussion 

We start with the zero temperature phase diagram in the plane U vs. 6n. An 
approximate analytic treatment is possible in this limit, which is used to check our 
numerical results. The eff'ect of temperature is considered afterwards. 

4.1. Zero temperature 

4.1.1. Approximate solution In Fig. [3] we represent the density of states of a biased 
bilayer with U = together with the low doping critical value Uc, as given by Eq. ( l27ll . 
In panel (b) of Fig. [3] a zoom in of the density of states close to the gap is shown. It 
is clear that the density of states diverges at the edge of the gap. As consequence the 
closer to edge of the gap the chemical potential is the lower will be the critical Uc value. 
This quantity is shown in panel (c) of Fig. [3] as function of the chemical potential fi and 
in panel (d) as a function of doping 6n. The lowest represented value of Uc is about 
Uc — 2.7 eV to which corresponds an electronic doping density 6n ~ 2.5 x 10~^ electrons 
per unit cell. The step like discontinuity shown in panels (c) and (d) for Uc occurs when 
the Fermi energy equals V/2, signaling the top of the Mexican hat dispersion relation. 
It is clear from panel (d) of Fig. [3] that in the low doping limit Uc is a linear function 
of doping 6n. This limit enables for an approximate analytic treatment which not only 
explains the linear behavior but also provides a validation test of our numerical results. 
Firstly we note that for very low doping the density of states in Eq. ( l27ll is close to the 
gap edge, |/i| ~ ^g/2, where A^, is the size of the gap Eq. ([32l) . In this energy region 
the density of states has a ID like divergence. [14] behaving as. 



VW -^s/2 

Using this approximate expression to compute the doping, 6n oc sign(/i) x Jj^ ,2 dx Pb{x), 
we immediately get 5n oc sign(/i)/pfe(/i) and thus Uc oc \5n\. In order to have an analytic 
expression for Uc in the low doping limit we have to take into account the proportionality 
coefficient in Eq. fl33l) . After some algebra it can be shown that the density of states 
per spin per lattice point near the gap edge can be written as. 



where x = [{^l + V'^)/{At^)Y/'^ , with F{x) and K(m) as in Eqs. ([HD and ^. The 
doping (5n, measured with respect to half filling in units of electrons per unit cell, can 
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Figure 3. (Color online) (a) - Density of states p(/i) per unit cell per spin of the 
bilayer problem with U = 0. (b) - Zoom of (a) near the gap region, (c) - Critical 
value Uc for ferromagnetism in the low doping, Sn, regime, (d) - The same as in (c) 
as a function of doping. The parameters are t = 2.7 eV, t± = 0.2t, V = 0.05 eV. The 
edge of the gap is located at Ag/{2t) ~ 0.00922. 



be written as, 



6n = sign(/i) x 8 / dx ph{x) 

'Ag/2 



tH 



-^/ Fix) ^IfwJv'^'"^^/'- 



(35) 



Inserting Eq. (l34ll into Eq. (l27l) . and taking into account Eq. (l35l) . we are able to write, 

-2 



Ur 



t TX 



Fix) 



. ^Fix)J 



6n . 



(36) 



A.(ti + v^) 

In panel (d) of Fig. [3] both the numerical result of Eq. f l27l l and the analytical result of 
Eq. (l36ll are shown. The agreement is excellent. 



4-1-2. Self- consistent solution We now need to solve the mean field equations in order 
to obtain the zero temperature phase diagram of the biased bilayer. In order to achieve 
this goal we study how m, Am, and An depend on the interaction U, for given values 
of the electronic doping 6n. 
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Figure 4. (Color online) Panels (a), (b), and (c) show the zero temperature self- 
consistent solution for m, Am, and An, respectively. The zero temperature phase 
diagram of the biased bilayer in the U vs. Sn plane is shown in panel (d). Symbols 
in panel (d) are inferred from panel (a) and signal a first-order phase transition; the 
solid [Eq. (|26l) ] and dashed [Eq. ((36|) ] lines stand for a second- order phase transition. 
The constant parameters are V = 0.05 eV, t± = 0.2t, and t = 2.7 eV. 



In panel (a) of Fig. [4] it is shown how m depends on U for different values oi6n. The 
chosen values of Sn correspond to the chemical potential being located at the divergence 
of the low energy density of states. The lower the Sn is the more close to the gap edge is 
the chemical potential and therefore the larger the density of states is. As a consequence, 
m presents a smaller critical Uc value for smaller Sn values. It is interesting to note that 
the magnetization saturation values correspond to full polarization of the doping charge 
density with m = Sn, also found within a one-band model. [46] In panel (b) of Fig. [4] we 
plot the Am mean field parameter. Interestingly the value of Am vanishes at the same 
Uc as m. For finite values of m we have Am > m, which means that the magnetization 
of the two layers is opposite. We therefore have two ferromagnetic planes that possess 
opposite and unequal magnetization. In panel (c) of Fig. [4] we show the value of An 
as function of U. It is clear that \Sn\ < \An\, which implies that the density of charge 
carriers is above the Dirac point in one plane and below it in the other plane. This 
means that the charge carriers are electron like in one plane and hole like in the other. 

In panel (d) of Fig. [4] we show the phase diagram of the system in the U vs. Sn plane. 
Symbols are inferred from the magnetization behavior in panel (a). They signal a first- 
order phase transition when m increases from zero to a finite value [see panel (a)]. The 
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Figure 5. (Color online) Effect of t±_ and V on the zero temperature Uc vs. 8n phase 
diagram: (a) fixed t\_ — 0.2i and varying V\ (b) fixed V = 0.05 eV and varying t^. 
For a given in the ferromagnetic phase establishes for U > Uc and the paramagnetic 
phase for U <Uc- 



full (red) line is the numerical self-consistent result of Eq. fl26l) . and the dashed (blue) 
line is the approximate analytic result given by Eq. ( l36ll . The discrepancy between lines 
and symbols has a clear meaning. In order to obtain both Eqs. ( l26ll and ( l27ll we assumed 
that a second-order phase transition would take place, i.e., the magnetization m would 
vanish continuously when some critical Uc is approached from above. This is not the 
case, and the system undergoes a first-order phase transition for smaller U values than 
those for the second-order phase transition case. There are clearly two different regimes 
in panel (d) of Fig. [H one at densities lower than 5n <lx 10~^, where the dependence 
of 5n on Uc is linear, and another regime for 5n > Ix 1Q~^ where a plateau like behavior 
develops. This plateau has the same physical origin as the step like discontinuity we 
have seen in panels (c) and (d) of Fig. [3l Clearly, as the density 5n grows the needed 
value of Uc for having a ferromagnetic ground state increases. This is a consequence of 
the diverging density of states close to the gap edge. As regards the limit 5n -^ Q \i is 
obvious from panel (d) of Fig.[4]that we have Uc -^ 0. It should be noted, however, that 
lowering the density 6n leads to a decrease of m and Am, as can be seen in panels (a) 
and (b) of Fig. [H Therefore, even though we have t/c — > in the limit (5n — > 0, we have 
also m — > and Am -^ 0, which implies a paramagnetic ground state for the undoped 
{6n = 0) biased bilayer. Only An remains finite at zero doping, in agreement with the 
observations that screening is still possible at the neutrality point {6n = ) . [49l [T31 150] 
So far we have analyzed the system for fixed values of the bias voltage, V, and 
interlayer coupling, t±. In Fig. [5] we show the effect of the variation of these two 
parameters on the zero temperature phase diagram. In panel (a) we have fixed the 
interlayer coupling, t± = 0.2t, and varied the bias voltage, ^(eV) = {0.01,0.05,0.1}; in 
panel (b) we did the opposite, with V = 0.05 eV and tj_/t = {0.05, 0.1, 0.2}. Essentially, 
raising either V or tj_ leads to a decrease of the critical interaction, Uc, needed to establish 
the ferromagnetic phase for a given Sn. The order of the transition, however, remains 
first-order: for a given 6n, the critical interaction Uc predicted by Eq. ( l26ll . which is 
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valid for a second-order phase transition, is always higher than what is obtained by 
solving self-consistently the mean-field equations, meaning that a first-order transitions 
is occurring at a lower Uc- It is interesting to note that the eff'ect of V and t± on the 
first-order critical Uc line is similar to what is expected for the usual Stoner criterion, 
where increasing either V or t± gives rise to an increase in the density of states at the 
Fermi level and a lower Uc thereof. 

The bias voltage and the interlayer coupling have also interesting effects on the 
magnetization, m, and spin polarization difference between layers, Am. Decreasing t± 
leads to a decrease in Am, and below some t± we have Am < m, as opposed to the 
case discussed above {V = 0.05 eV and t± = 0.2t). In particular, for V = 0.05 eV, we 
have already found Am < m for t± < O.lt. A similar effect has been observed when V 
is increased. For t± = 0.2t we have found Am < m for \^ > 0.1 eV. It should be noted, 
however, that m and Am are ^/-dependent. Increasing U leads m to saturate while Am 
keeps growing, as can be seen in panels (a) and (b) of Fig. [4] for the particular case of 

V = 0.05 eV and t± = 0.2t. This means that, depending on the value of the parameters 

V and t±, we can go from Am < m to Am > m just by increasing the interaction 
strength U. It can also be seen in panel (a) of Fig. [4] that m is completely saturated 
at the transition for 6n < dric ~ 6 x 10^^ electrons per unit cell, while for 5n > dric it 
saturates only at some U > Uc. Even though this behavior seems to be general for any 

V and t±, the value of Suc is not. In particular, we have found Suc to depend strongly 
on l^ - it seems to vary monotonically with V, increasing when V increases. Let us 
finally comment on the effect of V and t± on the charge imbalance between planes, An. 
Irrespective of V and t± we have always observed |5n| < |An|, which means that charge 
carriers are always electron like in one plane and hole like in the other. As expected, 
increasing/decreasing either V or t± leads to an increase/decrease of An. 

4.1.3. Understanding the asymmetry between planes The asymmetry between planes 
regarding both charge and spin polarization densities can be understood based on the 
Hartree-Fock bands shown in Fig. [6l The figure stands for V = 0.05 eV and tj_ = 0.2t, 
but can easily be generalized for other parameter values. 

It should be noted firstly that in the biased bilayer the weight of the wave 
functions in each layer for near-gap states is strongly dependent on their valence band 
or conduction band character. [TSl SH [50] Valence band states near the gap have their 
amplitude mostly localized on layer 2, due to the lower electrostatic potential —V/2 
[see Eq. ([4|)]. On the other hand, near-gap conduction band states have their highest 
amplitude on layer 1, due to the higher electrostatic potential -\-V/2 for this layer [see 
Eq. m- 

The case U < Uc shown in Fig. [6| (left) stands for the paramagnetic phase. The 
values m = and Am = seen in this phase are an immediate consequence of the 
degeneracy of up and down spin polarized bands. The presence of a finite gap, however, 
leads to the abovementioned asymmetry between near-gap valence and conduction 
states. As a consequence, a half-filled bilayer would have n2 = (4 + An)/2 electrons per 
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Figure 6. (Color online) Hartree-Fock bands for up (full lines) and down (dashed 
lines) spin polarizations. Three different cases are considered (from left to right): 
U <Uc,U> Uc, and U > Uc. 



unit cell on layer 2 (electron like charge carriers) and ni = (4 — An)/ 2 electrons per 
unit cell on layer 1 (hole like charge carriers), with An ^ 0. Even though the system 
studied here is not at half-filling, as long as \bn\ < \An\ the carriers on layers 1 and 2 
will still be hole and electron like, respectively. Note that as U is increased the charge 
imbalance An is suppressed in order to reduce the system Coulomb energy, as can be 
seen in panel (c) of Fig. [H From the band structure point of view a smaller An is the 
result of a smaller gap A^, which means that increasing U has the effect of lowering the 
gap. 

Let us now consider the case U > Uc shown in Fig. [6] (center). The degeneracy 
lifting of spin polarized bands gives rise to a finite magnetization, m 7^ 0. Interestingly 
enough, the degeneracy lifting is only appreciable for conduction bands, as long as 
U is not much higher than Uc- This explains why the total polarization m, and the 
difference in polarization between layers Am have similar values, m ~ Am, as shown in 
panels (a) and (b) of Fig. [4]- as only conduction bands are contributing to Am, the spin 
polarization density is almost completely localized in layer 1, where mi = (m+Am)/2 ^ 
m, while the spin polarization in layer 2 is negligible, m2 = (m — Am)/2 ^ 0. 

It is only when U ^ Uc that valence bands become non-degenerate, as seen in Fig. [6] 
(right). This implies that near-gap valence states with up and down spin polarization 
have different amplitudes in layer 2. As the valence band for down spin polarization 
has a lower energy the near-gap valence states with spin down have higher amplitude in 
layer 2 than their spin up counterparts. Consequently, the magnetization in layer 2 is 
effectively opposite to that in layer 1, i.e.. Am > m. This can be observed in panels (a) 
and (b) of Fig. [U where as U is increased the magnetization of the two layers becomes 
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Figure 7. (Color online) Panels (a), (b), and (c) show the finite temperature self- 
consistent solution for m, Am, and An, respectively, with temperature measured in K. 
The finite temperature phase diagram of the biased bilayer in the C/ vs. T plane is 
shown in panel (d). The constant parameters are V = 0.05 eV, t± = 0.2t, t = 2.7 eV, 
and 5n = 0.00005 e~/unit cell. 



opposite. 

We note, however, that the cases U > Uc and U ^ Uc mentioned above 
are parameter dependent. For instance, the valence bands can show an appreciable 
degeneracy lifting already for U > Uc, especially for small values of the t± parameter 
it_L < 0.05t). In this case the magnetization of the two layers is no longer opposite, with 
Am < m. This can be understood as due to the fact that as t± is decreased the weight 
of near-gap wave functions becomes more evenly distributed between layers, leading not 
only to a decrease in An but also in Am. As U is further increased the energy splitting 
between up and down spin polarized bands gets larger, enhancing Am. For U ^ U^ 
and depending on the parameters V and t_L, the magnetization of the two layers may 
become opposite even for small t^ values. 



4.2. Finite temperature 

Next we want to describe the phase diagram of the bilayer in the temperature vs. on-site 
Coulomb interaction U plane. This is done in Fig. [7] for a charge density (5ra = 5 x 10~^ 
electrons per unit cell. For temperatures ranging from zero to T =1.1 K we studied the 
dependence of m, Am and An on the Coulomb on-site interaction U . First we note that 
the minimum critical value Uc is not realized at zero temperature. There is a reentrant 
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behavior which is signaled by the smallest Uc for T = 0.06 ± 0.02 K. For temperatures 
above T ^ 0.1 K we have larger Uc values for the larger temperatures, as can be seen 
in panel (a). The same is true for Am, panel (b). As in the case of Fig. [H the value of 
Am, at a given temperature and U value, is larger than m. Also the value of An, shown 
in panel (c), is larger than 6n. Therefore we have the two planes presenting opposite 
magnetization and the charge carriers being hole like in one graphene plane and electron 
like in the other plane. In panel (d) of Fig. [7| we present the phase diagram in the T 
vs. U. Except at very low temperatures, there is a linear dependence of T on Uc- It is 
clear that at low temperatures, T ~ 0.2 K, the value of Uc is smaller than the estimated 
values of U for carbon compounds. [STl 152] 

4-3. Disorder 

Crucial prerequisite in order to find ferromagnetism is a high DOS at the Fermi energy. 
The presence of disorder will certainly cause a smoothing of the singularity in the DOS 
and the band gap renormalization, and can even lead to the closing of the gap. We 
note, however, that for small values of the disorder strength the DOS still shows an 
enhanced behavior at the band gap edges. [53l [54] The strong suppression of electrical 
noise in bilayer graphene [55] further suggests that in addition to a high crystal quality 
- leading to remarkably high mobilities[5^ - an effective screening of random potentials 
is at work. Disorder should thus not be a limiting factor in the predicted low density 
ferromagnetic state, as long as standard high quality BLG samples are concerned. 

Let us also comment on the next-nearest interlayer-coupling 73, which in the 
unbiased case breaks the spectrum into four pockets for low densities. [H] In the biased 
case, 73 still breaks the cylindrical symmetry, leading to the trigonal distortion of the 
bands, but the divergence in the density of states at the edges of the band gap is 
preserved. [54] Therefore, the addition of 73 to the model does not qualitatively change 
our result. 

5. Summary 

We have investigated the tendency of a biased bilayer graphene towards a ferromagnetic 
ground state. For this, we used a mean-field theory which allowed for a different carrier 
density and magnetization in the two layers. We have found that in the ferromagnetic 
phase the two layers have unequal magnetization and that the electronic density is 
hole like in one plane and electron like in the other. We have also found that at zero 
temperature, where the transition can be driven by doping, the phase transition between 
paramagnetic and ferromagnetic phases is first-order. 
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